home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_5.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.9 KB  |  185 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.5 (Newton-Raphson Iteration).
  9. % Section    2.4,    Newton-Raphson and Secant Methods, Page 84
  10. echo on; clc; format long; hold off; clear
  11. % This program implements the Newton-Raphson method.
  12.  
  13. %    Define and store the functions  f(x)  and  f'(x)
  14.  
  15. %    in the M-files  f.m  and  df.m   respectively.
  16. % function y = f(x)
  17. % y = x.^3 - x - 3;
  18.  
  19. % function y1 = df(x)
  20. % y1 = 3*x.^2 - 1;
  21.  
  22.  
  23. delete f.m
  24. diary f.m; disp('function y = f(x)');...
  25.            disp('y = x.^3 - x - 3;');...
  26. diary off;
  27.  
  28.  
  29. delete df.m
  30. diary df.m; disp('function y1 = df(x)');...
  31.             disp('y1 = 3*x.^2 - 1;');...
  32. diary off;
  33.  
  34. % Remark. f.m, df.m and newton.m are used for Algorithm 2.5
  35. f(0); df(0); % Test for files f.m, df.m
  36. pause    % Press any key to see the graph y = f(x).
  37.  
  38. clc; clg;
  39. a = -3.0;
  40. b = 3.0;
  41. c = -10;
  42. d = 10;
  43. h = (b-a)/150;
  44. X = a:h:b;
  45. Y = f(X);
  46. axis([a b c d]);...
  47. plot(X,Y,'-g');...
  48. hold on;...
  49. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  50. xlabel('x');...
  51. ylabel('y');...
  52. title('Graph of y = f(x).');...
  53. grid;...
  54. axis;...
  55. hold off;...
  56. shg; pause    % Press any key to perform Newton-Raphson iteration.
  57.  
  58. clc;
  59. %    Place the starting value in  p0
  60.  
  61. %    Place the abscissa tolerance in  delta
  62.  
  63. %    Place the ordinate tolerance in  epsilon
  64.  
  65. %    Place the maximum number of iterations in  max1
  66.  
  67. p0 = 2.0;
  68. delta = 1e-12;
  69. epsilon = 1e-12;
  70. max1  = 50;
  71.  
  72. [p0,y0,err,P] = newton('f','df',p0,delta,epsilon,max1);
  73.  
  74. pause    % Press any key for the list of iterations.
  75.  
  76. clc; clg;
  77. a = 0.9;
  78. b = 2.1;
  79. h = (b-a)/150;
  80. X = a:h:b;
  81. Y = f(X);
  82.  
  83. max1 = length(P);
  84. clear Vx Vy
  85. for i = 1:max1,
  86.   k1 = 2*i-1;
  87.   k2 = 2*i;
  88.   Vx(k1) = P(i);
  89.   Vy(k1) = 0;
  90.   Vx(k2) = P(i);
  91.   Vy(k2) = f(P(i));
  92. end
  93. a = 1.6;
  94. b = 2.05;
  95. c = -0.5;
  96. d = 3.5;
  97. Z0 = zeros(1,length(P));
  98. axis([a b c d]);...
  99. plot(X,Y,'-g',Vx,Vy,'-r',P,Z0,'or');...
  100. hold on;...
  101. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  102. xlabel('x');...
  103. ylabel('y');...
  104. title('Graphical analysis for Newton`s method.');...
  105. grid;...
  106. axis;...
  107. hold off;...
  108. shg; pause    % Press any key to continue.
  109.  
  110.  
  111. J = 1:max1;
  112. Yp = f(P);
  113. points = [J;P';Yp'];
  114. Mx1 = 'Iterations for Newton`s method.';
  115. Mx2 = '     k                  p(k)               f(p(k))';
  116. Mx3 = 'The solution is:';
  117. Mx4 = 'The error estimate for p is  ± ';
  118. clc,echo off,diary output,...
  119. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  120. disp('Iteration converged quadratically to the root.'),...
  121. disp(''),disp(Mx3),disp(''),disp('p = '),...
  122. disp(p0),disp('f(p) = '),disp(y0),...
  123. disp([Mx4,num2str(err)]),diary off,echo on
  124.  
  125. pause    % Press any key to perform Newton-Raphson iteration.
  126.  
  127. clc;
  128. %    Place the starting value in  p0
  129.  
  130. %    Place the abscissa tolerance in  delta
  131.  
  132. %    Place the ordinate tolerance in  epsilon
  133.  
  134. %    Place the maximum number of iterations in  max1
  135.  
  136. p0 = 0.0;
  137. delta = 1e-12;
  138. epsilon = 1e-12;
  139. max1  = 12;
  140.  
  141. [p0,y0,err,P] = newton('f','df',p0,delta,epsilon,max1);
  142.  
  143. pause    % Press any key for the list of iterations.
  144.  
  145. a = -3.5;
  146. b =  0.5;
  147. c = -30;
  148. d = 5;
  149. h = (b-a)/100;
  150. X = a:h:b;
  151. Y = f(X);
  152. max1 = length(P);
  153. clear Vx Vy
  154. for i = 1:max1,
  155.   k1 = 2*i-1;
  156.   k2 = 2*i;
  157.   Vx(k1) = P(i);
  158.   Vy(k1) = 0;
  159.   Vx(k2) = P(i);
  160.   Vy(k2) = f(P(i));
  161. end
  162. Z0 = zeros(1,length(P));
  163. axis([a b c d]);...
  164. plot(X,Y,'-g',Vx,Vy,'-r',P,Z0,'or');...
  165. hold on;...
  166. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  167. xlabel('x');...
  168. ylabel('y');...
  169. title('Graphical analysis for Newton`s method.');...
  170. grid;...
  171. axis;...
  172. hold off;...
  173. shg; pause    % Press any key to continue.
  174.  
  175.  
  176. J = 1:max1;
  177. Yp = f(P);
  178. points = [J;P';Yp'];
  179. Mx1 = 'Iterations for Newton`s method.';
  180. Mx2 = '     k                  p(k)               f(p(k))';
  181. clc,echo off, diary output,...
  182. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  183. disp('Iteration did not occur.  This is a case of "cycling."'),...
  184. diary off, echo on
  185.